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Abstract 

The distribution of visible matter in the universe, such as galaxies and galaxy clusters, has 
its origin in the week fluctuations of density that existed at the epoch of recombination. The 
hierarchical distribution of the universe, with its galaxies, clusters and super-clusters of galaxies 
indicates the absence of a natural length scale. In the Newtonian formulation, numerical simulations 
of a one-dimensional system permit us to precisely follow the evolution of an ensemble of particles 
starting with an initial perturbation in the Hubble flow. The limitation of the investigation to 
one dimension removes the necessity to make approximations in calculating the gravitational field 
and, on the whole, the system dynamics. It is then possible to accurately follow the trajectories of 
particles for a long time. The simulations show the emergence of a self-similar hierarchical structure 
in both the phase space and the configuration space and invites the implementation of a multifractal 
analysis. Here, after showing that symmetry considerations leads to the construction of a family of 
equations of motion of the one-dimensional gravitational system, we apply four different methods 
for computing generalized dimensions Dq of the distribution of particles in configuration space. We 
first employ the conventional box counting and correlation integral methods based on partitions of 
equal size and then the less familiar nearest-neighbor and k-neighbor methods based on partitions 
of equal mass. We show that the latter are superior for computing generalized dimensions for 
indices q < —1 which characterize regions of low density. 

* b.miller@tcu.edu 
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I. INTRODUCTION 


The analysis of large scale surveys shows that the grouping of visible matter presents a 
hierarchical structure on very large scales [1]. Stars are grouped into galaxies, galaxies into 
clusters, clusters into super-clusters, and so on. The shape of the ultimate element of this 
hierarchy has yet to be determined. Recently, a new grouping technic using the peculiar 
velocities of galaxies led to the discovery of the Laniakea supercluser, thereby extending the 
structure horizon [2]. This hierarchical regrouping, for which the enlargement of a structure 
reveals the existence of a smaller one, suggests a fractal arrangement of matter[3, 4]. A true 
mathematical fractal object exhibits structures on all scales [5]; however, in practice, for the 
case of a physical object, there are both upper and lower bounds to the scaling structure 
[6]. A quantitative measure of the range of scales that accommodate a hierarchical (fractal) 
distribution in the universe is provided by the correlation function of the distribution of 
galaxies obtained from recent surveys [3]. 

Numerical simulation has made a significant contribution to the study of structure for¬ 
mation. It allows us to follow the dynamical evolution of these structures that change 
too slowly to be observed [7]. A number of specihc, three dimensional, hydrodynamic and 
N-body codes are employed. However they require considerable resources and are limited 
in their ability to reveal fractal structure by the hnite resolution that can be realized [8]. 
Alternatively, by limiting simulations in jj, space to one space dimension and an A^-body 
description, we can increase the number of particles per dimension and treat the dynamics 
exactly, thereby retaining all the information concerning the particle trajectories. Of course, 
the hnite size of the memory that represents numbers on computers doesn’t allow the strict 
reversibility of the trajectories. 

A purely one dimensional model dates from the year 1990 [9]. The spherical version 
(hereinafter the Q model) was introduced by Fanelli and Aurell [10]. An inhnite version of 
the system was studied by Gabrielli et al. including the temporal evolution of the power 
spectra [11]. The multifractal properties of a periodic version have been studied by Miller 
and Rouet [12, 13]. Joyce et al. have studied the virialization of the clusters with regard 
to three dimensional observations. Virialization of gravitational systems have been given 
in numerous studies, for example [14, 15]. The fractal properties of these structures have 
notably been placed in evidence by Tsuchiya et al. [16]. 
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In the earliest studies of the one-dimensional system [9], the initial condition was chosen 
as a “water bag” in which the particles were equally spaced in position, but their velocities 
were chosen independently and at random from a uniform distribution symmetric about the 
origin. After the system evolved, a fractal analysis was performed using the box-counting 
dimension. Depending on the Jeans’ length associated with the initial condition, hierarchical 
scaling was revealed for the mass distribution in both phase space and conhguration space. 
Our more recent collaborative efforts have followed three principle avenues: the investigation 
of the possible existence of multifractal structures, the influence of scale free initial condi¬ 
tions closer to those following inflation revealed by WMAP, and the influence of changes 
in the parameters of the one-dimensional model [12, 17]. In so doing we also showed how 
to rigorously formulate the evolution of a one-dimensional self-gravitating system obeying 
periodic boundary conditions [13]. 

In the present work, after considering how symmetry dictates the construction of the 
equation of motion of the one-dimensional gravitational system, we focus on a comparison 
of different methods for carrying out the fractal analysis of the resulting distribution in 
conhguration space. A development of our recent work is that the standard methods for 
computing generalized, or Renyi, dimensions Dg exhibit problems for the low density regions 
characterized by negative values of q [12]. Here we hrst present a derivation of the class of 
one-dimensional models that can be constructed based on symmetry arguments. We then 
select the model originally introduced by Rouet and Feix which is the most self-consistent 
of the class of one-dimensional models. After following the evolution from small initial 
perturbations in the Hubble how to a highly clustered state, we compare mass oriented 
methods of fractal analysis with the results of more standard approaches based on partitions 
of equal size. Following this introduction, in section H, the family of models and algorithms 
for following the dynamics of the particles which compose the system are presented. A typical 
simulation is presented in section HI along with the results of the multifractal analysis. These 
results are discussed in the last section of the article. 


II. MODELS AND ALGORITHM 

In the following we consider a segment of the universe with extension sufficiently small 
for the Newtonian approximation to remain valid. The system will be described by a set of 
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N particles interacting pair-wise according to Newton’s law of gravity. We assume that the 
distribution of matter is highly symmetric. Specihcally, in three dimensional space, three 
models are considered for which the distribution symmetry is either planar, cylindrical or 
spherical. Let x refer to the coordinate of a particle {x corresponds to the radius of a cylinder 
in the cylindrical case for example). Its equation of motion is written 


d^x 

dt^ 


E{x, t) 


( 1 ) 


where E{x,t) is the gravitational field. To accommodate the expansion two new variables 
of space and time are introduced [18] where C(t) represents the cosmological scale factor [1] 


X = C(t)x 

dt = A‘^{t)di (2) 


where 


C{t) = 

A^{t) = (3) 


wjq is dehned as the Jeans frequency, where po is the mass density taken at 

the initial time f = to at which the size of the physical and rescaled universe correspond. So 
C(to) = ^(to) = 1 and thus 


to — 


l^Jo 


( 4 ) 


Assuming that the held in the new variables follows the usual Poisson law, the equation 
of motion in the rescaled variables is written: 


d^x 

dt^ 


+ 2A^ 






( 5 ) 


where C notes the derivative of C with respect to time t and C the second derivative. In 
taking a = 2/3 and /3 = 1, equation (5) no longer contains any time dependent coefficients. 
It becomes 


d^i; 

dt^ 


1 da; 

+ - 


2 

9 


(7^ajo)^£ 


E[x^ i) 


( 6 ) 
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With the choice a = 2/3, C{t) represents the Hubble expansion for a matter dominated 
universe that is appropriate for the epoch just following recombination and this is what we 
assume here [1]. Then the new variables represent a comoving frame and only the residual 
movement will be calculated. Note that the new unit of time dt is now constant: 


uj(t)dt = a/ 4:nGp(t)dt (7) 

= ( 8 ) 

The transformation (2) reveals two new terms in the equation (5), a friction and a force 
proportional to x. In the following, it will be shown that the coefficient 7 depends on the 
type of symmetry and that, with one exception, the latter term of the left hand side of 
equation (5) must be modihed so as to neutralize a uniform system that exactly follows 
the expansion. Under these conditions, the rescaled system will be static. It is a planar 
perturbation of this static universe that we follow hereafter. 

A. Symmetries and dimensions 

From the imposed symmetry, i.e., given that the field depends only on x, in d dimensions 
the divergence operator is written: 

1 dx'^~^E 

where d = 1 for planar symmetry, d = 2 for cylindrical symmetry, d = 3 for spherical 
symmetry, etc.. From the Poisson equation, in the comoving frame the scaled field then 
becomes 


4 

E = —-TiGpx. (10) 

for the mean Hubble flow as p is constant. In that case the hrst two terms of equation ( 6 ) 
vanish. By requiring that the force is proportional to x, since the transformation compen¬ 
sates for that of gravity, we have 


-^'^Gpx = --{juj^fx 

and therefore 


( 11 ) 
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p = Po and 7 = 


v^' 


( 12 ) 


Equation (6) can now be written in arbitrary dimension d as 


d^x 1 da; 1 o - 


(13) 


d«i ° di <i *' 

Considering our earlier assumption that the static state (mean flow) is excited by a planar 
perturbation of amplitude y we have 


= £(;),£) (14) 

The held is now that of a planar system so, in order to satisfy global neutrality (background 
+ held), the coefficient 1/d has been removed from the background term. When d = 1, this 
seemingly arbitrary change is not necessary because both the disturbance and the system are 
planar. Following the notation of Fanelli [10] this model (d = 1) will be called the RF model 
and the spherical model (d = 3 ) the Q model (or Quintic). Note that, for a spherically 
symmetric system of dimension d, the relation (9) remains valid. If d approaches infinity, 
the coefficient of friction approaches 0. While this requires imagining a universe of more 
than 3 space dimensions, this observation gives a physical meaning to the model without 
friction, or H-model, which has also been studied [12] elsewhere. 

The system is composed of N particles which are infinite parallel planes themselves, 
with equal and constant surface mass density p, (with the effect of expansion taken into 
account by the scaling). They experience friction and are bathed in a fixed, neutralizing, 
homogeneous background. In this sense, this system is the inverse of a plasma, since here the 
particles are attracted to each other and repulsed by the background medium. In cosmology 
it is customary to consider a segment of the universe that is small compared to the Hubble 
distance, but large enough to contain many galaxies and, optimistically, the largest observed 
structures. Since the size is much less than the Hubble distance, relativistic effects can be 
ignored and Newtonian mechanics suffices to follow the evolution [19]. In order to minimize 
edge effects and more closely represent a segment of an infinite system, periodic boundary 
conditions (PBC) are assumed to apply [20]. Analogously, for our one dimensional models, 
we consider a slice of the universe that also obeys periodic boundary conditions. Then the 
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gravitational field due to a single plane must also include the contribution from the sum of 
all of its inhnite replicas, commonly referred to in the literature as an Ewald sum [13]. 

Since the held due to an isolated sheet of mass is constant in space, obtaining the correct 
held for PBC is not a trivial problem. Elsewhere we have shown that each mass sheet (or 
particle) carries with it its own neutralizing background [13]. This turns out to be the same 
background contribution that appears in the transformation to the comoving frame for the 
RE model without further adjustment. In this sense, PBC are the only correct boundary 
conditions that are compatible. The held experienced by a particle in the primitive cell 
satisfying PBC is given by 


E{x) = inmG 



+ ^(Nk{x) - iVi(x)) 


(15) 


where 2L is the system size and Xc is the system center of mass computed in the primitive cell 
[13]. The role of Xc is to insure that no interruption of the held is experienced by an interior 
particle when a diherent particle crosses a cell boundary and enters from the opposite side. 

Operationally, one way this Ewald sum can be achieved is by polarizing the system 
boundaries: each side acts as a reservoir of initially neutral particles that can be loaded 
with particles of ehectively “positive” and “negative” mass. As a particle reaches an edge, 
this edge becomes positively “charged” and the other one is negatively “charged”, thus 
introducing a compensating dipolar gravitational held. This technique is usually used in 
the case of a one-dimensional plasma and insures that there is no discontinuity in the held 
when a particle crosses a cell boundary. As explained in [13] it is thus possible to follow 
the particle trajectories exactly until they intersect. Then, rather than the usual pattern 
of dynamical evolution, which is to advance molecular particles along their trajectory ac¬ 
cording a hxed time step, here we adopt an event driven scheme where the particles are 
advanced to their next crossing time with their neighbor. As the solutions of the equations 
of motion are known, it is possible to calculate the shortest crossing time of a particle with 
its neighbor and then restart the process. For the RE (Quintic) case, this means solving 
a cubic (hfth order) equation and determining the root that yields the smallest positive 
crossing time. Thus the RE model has the advantage that the crossing times can be found 
analytically. For diagnostics, the particles are temporarily advanced to the current time. 
This is straightforward as the analytic solutions of the equations of motion between crossing 
events are known exactly. Typically the trajectories of A^ = 2^® particles are followed in 
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quadruple precision. This precision is necessary to follow a large number of particles for a 
sufficient time to enable large size structures to grow, while retaining a high precision for 
the distribution of structures of smaller size (cf. hgure 1). 


III. MULTIFRACTAL ANALYSIS OF SIMULATIONS 

Following inflation, density fluctuations in the universe can be modeled as a Gaussian 
random field with a scale-free (power law) power spectrum [21], In the three dimensional 
universe, the exponent corresponding to a scale-free potential is unity [21], Initial conditions 
for 3D simulations of the expanding universe are guided by these considerations. For the 
simulation data reported here, we selected the RF model and, in order to achieve scale-free 
potential fluctuations in ID, the power spectrum of the density fluctuations at the initial 
time is chosen to vary as k^, where k is the wave number. The construction of the initial 
distribution of the particles allowing us to sample this spectrum is given by Miller et ah 
[12]. In the past, other initial conditions were adopted, particularly the water bag model 
for which the particle velocities are distributed randomly, following a uniform distribution 
between ±a, while the positions are equally spaced. If the friction term is omitted, it is 
possible to write a dispersion equation which shows that the system is unstable for wave 
numbers k such that k > kc = ‘I'njXj where Aj is the Jeans length. This particular initial 
condition was used to set the Jeans length ‘I'najuj. Nevertheless, for the initially cold system 
simulated here, the Jeans length tends to 0, and these initial conditions, with a power law 
spectrum of the density, do not favor any particular length scale. However, care must be 
taken in the choice of the spectral index [21]. If we take k~^ for example, we will give too 
much energy to large scales and we will just arrive at one large cluster, perhaps with a few 
small ones, while we require /c < 4 to guarantee momentum conservation [21]. 

Figure 1 shows a sequence of snapshots of the density and distribution of iV ~ 2^^ particles 
in the phase space. The initial fluctuations cause the system to rapidly enter a non-linear 
regime wherein structures of small size appear. These structures locally regroup to form 
larger structures that are themselves regrouped, and so on. Two successive enlargements of 
a larger structure formed at time T = 20 are shown in hgure 2. It can be seen both in the 
phase space, but also on the density peaks, that these large structures are formed of smaller 
structures, themselves resulting from an assembly of structures of less size. 
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FIG. 1. Snapshots of the distribution of particles in configuration and phase space for a system of 
N = 8191 particles at the instants T = 0, 8, 14, 18 and 20. 

This hierarchical structure suggests that, at least over a finite range of scales, the system 
is fractal. To test this hypothesis, a multifractal analysis on the distribution of particles in 
the conhguration space is conducted using four methods. Fractals can be characterized by 
their dimension. For simple fractals, such as a Koch curve or the Cantor set, a single number 
is sufficient, whereas for systems that are less homogenous, a continuum of dimensions is 
required. Considering the large variation in density exhibited on the line and in the phase 
space in the results of our simulations, it would seem that the latter approach has to be 
pursued. 

Renyi introduced the idea of generalized dimension Dq by hrst partitioning the embedding 
space into Nb equal size cells and determining the measure /ij associated with each. For the 
case of our simulations we consider the local density fii = Ni/N within each box where Ni 
is the number of particles in this box. Let C{q,l) be the partition function identihed with 


9 




































































T=7() 





z.oom 1 


s 


4 - 


-2000 


250 



-170 
-1300 -2000 

224 



-1300 


1810 -1970 



1810 


FIG. 2. Two successive zooms of a patch of phase space (upper figure) at T = 20 . The enlarged 
areas are bounded by dashed lines. 


a particular decomposition: 


NiU) 

C{q,l) = /i-(0 and limC'(g,/) = P. 

i=l 

Then the generalized or Renyi dimensions are defined by 


Dq = lim 


Dn = lim 

^ Z-5-0 


1 ln(C(g,/)) 
1^0 q — 1 ln(/) 


HI) 


for g 7 ^ 1 
for g = 1 


(16) 


(17) 


The exponent Tq is related to the Renyi dimension of order q hj Dg = Note that, from 
the dehnition, regions of high density contribute strongly for positive values of g, and low 
density regions more strongly support negative values. 

In practice, for simulated or experimentally observed data, it is not possible to go to the 
limit that the dehnition demands. To compensate, two general methods for determining 
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fractal dimensions have evolved. In one, the space in which the fractal set is embedded is 
partitioned into subsets of equal size, whereas, in the second, it is partitioned into sets of 
equal measure or mass. The first method follows more closely from the dehnition of Renyi 
Dimension. It is realized in both the box counting and correlation methods, which are the 
most popular for estimating dimensions of natural sets. The second approach is realized in 
the near-neighbor [22] and k-neighbor [23] methods and has the advantage of only including 
occupied cells in the partition. Recently, to gain insight concerning their useful regimes, we 
have applied these methods to standard sets with well-characterized fractal properties [24]. 
To paint as clear a picture as is possible regarding the fractal nature of the gravitational 
simulations, here we will pursue each approach and see what further insights they provide. 

In the box counting (BC) method, since the limit in Eq. (17) cannot be performed in 
practice, rather plots of In(Cq) vs ln(/) are studied to determine if linearity is present, 
suggesting power law behavior in the partition function that is required for a well dehned 
limit. We operationally define the generalized dimensions using 


r 1 ln(C(g,/)) 

I O' — 1 ln(/) 

[- —ta(i)— 


q^l 
q = l 


(18) 


The exponent Tq is related to the Renyi dimension of order qhy Dq= The representation 
of hi{C{q,l))/{q — 1) as a function of ln(/) gives the value of Dq in the zone for which the 
plot is linear. Note that if / is too large, then the slope is equal to 1 which indicates that 
the system is homogeneous on this scale. Conversely if I is too small, then the slope is 
zero, which is an effect of the discretization of the system. Due to the simplicity of the 
method, it is widely used among researchers. However, it has been pointed out by many 
that this method and, more generally, methods that involve partitions of the same size such 
as the correlation method, do not work well for q < 1. [24, 25] In the correlation integral (Cl) 
method introduced by Grassberger and Procaccia [26], one fixes on a selected set of reference 
points and examines the distribution of points in neighborhoods of equal size surrounding 
each of them. This provides an alternative formulation for a partition function Iq{l) for a 
partition of the system with cells of equal size. Let 

1 ^ 
i=l 
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where, here, 


1 ^ 
i=i 


where 9{x) is the unit step function and 


Dq = 

^ g- 5 >o q — 1 ln(/) 


( 20 ) 


( 21 ) 


as before. Then the generalized dimension Dq can be obtained in similar fashion to the BC 
method. 

To formulate approaches based on partitions of equal measure or mass, let 6j{k,n) be 
the distance between a point xj chosen from our set and the nearest neighbor to Xj. 
Next construct the sum of the moments Sj{k,n) from a set of n reference points chosen at 
random: 


AW(A;,n) =-V57(A;,n). (22) 

n •' 
i=i 

Van der Water and Schramm have shown that [23], in the limit of large n, 

where o; is a constant independent of 7 and the Dimension Function, D{'j), can be thought 
of as alternative generalized dimension. It is related to the Renyi dimension through 

D[j={l-q)Dq]=Dq. (24) 

Note that the average of 6j from a single set is used in Eq. (22) whereas the derivation of 
Eq. (23) is based on the ensemble probability. 

Using their result, two different approaches for determining £*( 7 ), and therefore Dq, can 
be extracted. For a fixed value of k, the partition elements all consist of intervals with the 
same number of particles. One can then investigate the dependence of on the sample 
size n and obtain —1/D{'y) from the slope of plots of log(AO)/iF]) vs log(n). In doing so 
note that care must be taken to avoid the singularities of the Gamma function [27]. For 
the special case where k = 1 this is known as the nearest-neighbor (NN) method and was 
extensively investigated by Badii and Politti [22], later by Broggi [28] and, very recently, 
by us [24]. However, it has also been applied with hxed values of k as large as 300 [28]. 


0:^(7) 


F(fc + 7/D(7)) 
F(A;) 


1/7 


= k) 


(23) 
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Alternatively, we can fix the number of sample points n and investigate the scaling of A 
with k. In so doing we are considering a cumulative sequence of partitions based on elements 
of increasing mass or measure. This approach is known in the literature as the k-neighbor 
(KN) method. Because it incorporates a wide range of k values, it results in a more global 
measure of dimension and is less susceptible to local fluctuations that influence the other 
methods[24]. For large k, a simple approximate relation can be obtained: [23] 

[A(^)(A;, n)] ^ 7 ) (25) 

where G{k,'y) is a correction function and is close to unity for large k. The Dimension 
Function D{'y) can be estimated by setting G{k, 7 ) = 1 in the first iteration. The dependence 
of the correction function G{k, 7 ) on k and 7 can be obtained from Eq. (23) with the value 
of from the first iteration. The Dimension Function £*( 7 ) is then updated using this 
G{k,'-f). After a few iterations, the numerical results for £*( 7 ) typically converge to a single 
value for each 7 . 

The set of analyses are carried out at time T = 20 on a simulation with N = 2^® particles. 
The particle distribution is similar to that of Figures 1 and 2 for which the number of particles 
is 2^^, except that the number aggregate is now 32 times greater, improving the statistics. 
Regarding the BC method. Figure 3 shows a typical plot of ln(Cq)/(g — 1) as a function of 
ln(/) for different values of q. You can notice the intervals where Dg = 0 (zero slope) for 
I —)■ 0 and where Dg = 1 for large values of 1. Between these two zones, there is a third 
scaling interval which seems to expand with the value of q. For q = 0, for example, it varies 
from about 10“^ to 10^, encompassing six decades. The values of the slope for these scaling 
regions are identihed for several values of q. The function of Dq vs q thus obtained is shown 
in Figure 4. Similar results obtained by the Cl method are also shown in Figure 4. The two 
curves are similar. For g > — 1, they are decreasing, but increasing for g < — 1. However, 
while theory tells us that a general property of Dq is that it is not strictly increasing, this is 
not the case here. This situation led us to employ mass oriented partitions to have better 
insight, especially for negative values of g. 

As the nearest neighbor method is sensitive to the local statistical noise, we used k = 3 
instead of k = 1 , which we refer to as the “near-neighbor” method. As shown in Fig. 5, we 
have successful scaling for the positive range of 7 . However, as 7 increases, the contributions 
from a few sample points start to dominate the sum in Eq. 22, and therefore the scaling 
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ln(l) 


FIG. 3. BC method: ln{Cq)/{q — 1) vs. ln(Z) for several values of q. The slope of the linear central 
zone gives Dg (cf. hgure 4). 

becomes less stable. By increasing the fixed value of /c, the location of the singularity is 
shifted further into the negative range of 7 , thus increasing the size of the reliable range for 
this method. 

The evolution of generalized dimension Dg over time T computed with the near-neighbor 
method is shown in Fig. 6 . At the beginning of the simulation, we obtain a constant Dg 
with a value close to 1 as expected. For the range q > 2, the generalized dimension starts 
to diverge from the expected spectrum. Accordingly, we conclude that the near-neighbor 
method with k = 3 is not reliable for q > 2. This conforms with the analysis performed 
on sets with known fractal properties [24]. For the range where this method is reliable, we 
obtained a smooth, non-increasing, spectra for Dg. For the range q > —2.5, Dg generally 
increases over time whereas for q < —2.5, the spectrum decrease over time. It is noteworthy 
that the generalized dimension is almost invariant over the entire duration of the simulation 
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q 


FIG. 4. Generalized dimension Dg{q) obtained at T = 20 by the BC (squares) and IC (circles) 
methods for a simulation of = 2^® particles. 

at g = —2.5 with Dq = 1. This observation may not be accidental and warrants further 
study. 

To evaluate D{'j) using the k-neighbor method, A'^, which has the dimensions of length, 
is plotted vs. the number of neighbors k for several values of 7 on a log-log scale. From the 
slope, the values of Dq for differing q can be deduced following the relation (24). In Fig. 7, 
we show the plots obtained with the fc-neighbor method. We divide the range of k into four 
regions depending on the behavior of the weighted average of the neighbor distance. In 
this particular plot, we see large gaps in the plots for positive 7 at, for example, k = 10. The 
range fc < 10 is known to have singularities for the negative range of 7 , and therefore we 
typically do not include this region in the subsequent computation. The range 10 < /c < 100 
exhibits typical fine structures which are the hallmark of a fractal. The gaps correspond to 
lacunarity in the simulated set. The generalized dimension Dq extracted from this region is 
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FIG. 5. This plot shows the scalings with the near-neighbor method. The slope of the best-fit 
line for each 7 is taken as the Dimension Function, 


included in Fig. 8. Although the sample points are limited, and therefore a careful analysis is 
required, the spectrum is more mono-fractal like with Dg being significantly smaller than 1. 
The range 100 < k < 1200 clearly has different slopes for different values of 7 . The spectrum 
from this range is said to be multifractal as shown in Fig. 8. For the range k > 1200, all 
the plots with various 7 merge into a single line whose slope is close to 1. This shows that 
the simulated set is homogenous on large scales. In Fig. 8, we computed the generalized 
dimension using several different scaling ranges of k. We can see that the result is sensitive 
to the choice of scale. The result from the near-neighbor method is included for comparison. 
Here, we used an identical range of k for different 7 . While the spectra computed with some 
ranges more nearly resemble the spectrum with the near-neighbor method, we were not able 
to find a range of k where the generalized dimension closely traces the spectrum from the 
nearest neighbor method. For q> —1, the resulting values are compatible with those found 
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q 

FIG. 6. This plots shows the generalized dimension using the third nearest neighbor method for 
different times T. 

using the BC and IC methods. 

It is useful to check values of Dq by another method. For a fractal, the autocorrelation 
function of the density decays as a power law. From the autocorrelation function, or its 
Fourier transform which yields the power spectrum of the density fluctuations, it is possible 
to obtain the correlation dimension D 2 . We have the relation D 2 = —n where n is the 
power of the wave number k in the power spectrum of the density fluctuations [12]. The 
power spectrum is shown in Figure 12 at different times. At T = 0, the initial condition 
imposes a slope of 3 for the longer wavelengths, which is verihed in the hgure. For large 
wave numbers {k > kc = 27 r/l) the spectrum is 1, which reflects the lack of correlation 
between particles at separations of less than unity. Over time, one sees the formation of an 
intermediate scaling range that exhibits power law behavior. The size of the scaling region 
grows through decreasing its lower boundary and increasing the upper boundary: the slope 
is always the same. At T = 20, a linear regression analysis in the range .01 < fc < 1000 
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FIG. 7. This log-log plot shows the scalings with the A:-neighbor method at T = 20. 

yields a slope n = —0,45 and therefore D 2 = 0,45. This value is quite consistent with those 
found by the IC, BC and KN methods. 


IV. DISCUSSION AND CONCLUSION 

A simple one-dimensional model that only includes expansion and gravitation (Newto¬ 
nian) and is subject to periodic boundary conditions was investigated numerically. This 
model belongs to a family of one-dimensional models of expansion which only differ by their 
symmetry in the original space. The use of the quasi-invariance group [29], for which the 
law of transformation of the spatial variable is different from that of the time unit, provides 
an autonomous equation of motion. Other choices are possible, such as taking the same 
transformation law {A{t) = C{t)). We then obtain a more standard form of the comoving 
coordinates. Friction disappears, but then a gravitational constant which depends explicitly 
on the time is introduced. [1]. This family of models, which describe the evolution of a 
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FIG. 8. This plot the sensitivity of the generalized dimension Dq to the choice of the range from 
which the best-fit line is extrapolated in the fc-neighbor method. The selected ranges are shown 
in the legend. The generalized dimension Dq computed with the near-neighbor method is also 
inclnded for comparison. 


planar disturbance, differ only in the value of the coefficient of friction. This coefficient 
approaches zero if we consider a high-dimensional d space. Its properties were investigated 
elsewhere where it was referred to as the Hamiltonian case [12]. 

For three-dimensional simulations, it is necessary to introduce approximations in the 
gravitational held at both short and long distances. In particular, a short-range cut-off is 
necessary to maintain a manageable execution time. The inhuence of these approximations 
of the force law on the predicted fractal properties is not known, but could be signihcant. It is 
likely that they would limit the useful scaling range. In one dimension, these approximations 
are not necessary. The equations of motion are exactly integrated between particle crossings 
so it is possible to write an event driven algorithm to follow the motion with the precision 
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FIG. 9. This plot shows the generalized dimension Dq for selected times T, computed with the 
/c-neighbor method. 

imposed by the word-length maintained by the computer for a set of N particles. This is 
important for fractal analysis that requires us to keep all possible precision for the particle 
trajectories. Also, in contrast with higher dimension, elsewhere we have shown that it is 
possible to analytically evaluate Ewald sums in ID and thus exactly model periodic boundary 
conditions [13]. As we have shown, these are the only boundary conditions that are consistent 
with a comoving frame. 

The simulations were performed on an initially cold system with a density spectrum that 
follows a power law. This scale-free formulation is commonly used by astrophysicists [30] 
and is motivated by the generally accepted characteristics of the inflationary density field. 
Starting with a nearly homogeneous distribution of N particles, simulations show a successive 
fractionation and recombination resulting in the formation of a hierarchical structure of 
clusters in both phase space and on the real line. After a sufficiently long time, only a single 
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FIG. 10. This plot shows the values obtained for the generalized dimensions using different, 
mass-oriented, numerical methods. The corresponding value of 7 is also shown. 


structure will remain. This is only due to the finite number of particles. However, even if 
the boundary effects become important by then, the structure of the virialized cluster thus 
formed could be explored, pointing to differences with virialization of an isolated particle 
system which arises only as a consequence of gravity (without expansion, and therefore 
without friction or a continuous background). 

In earlier work we showed that the analysis by the BC and IC methods give a generalized 
dimension Dq{q) which increases when q increases for g < — 1, then remains almost constant 
for g > 1 [12]. However, a general property of the Dq curve is that it should be non-increasing 
[6]. Therefore, the results for g < —1 obtained by the BC and IC methods are questionable. 
As we have shown here, in contrast with the former methods based on partitions of equal size, 
the analysis based on partitions of equal mass or measure provides a generalized dimension 
Dq{q) which decreases as g in the low density regions. The fact that the results from the 
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FIG. 11. Generalized dimensions Dq{q) obtained at T = 20 by the fe-neighbor method (triangles). 
The results obtained by the BC and IG methods given in figure 4 are shown for comparison. 


two distinct mass-oriented methods do not show complete agreement may be explained 
by the hnite degree of hierarchy of the simulated set. Since the nearest neighbor method 
(NN) is concerned only with the local statistics of the given set, the method may extract 
generalized dimensions associated with its subsets that can have different properties from 
its global properties. From the scaling plots of the /c-neighbor method, we can see that 
scaling is observed in limited ranges which do not generally include small values, so k = 3 
is avoided. This may also contribute the the difference between the KN and near-neighbor 
methods since the latter explicitly employed k = 3. For the high density regions, three 
methods, namely, the BC, IC and near-neighbor, are known to work well and they agree for 
g > — 1. In addition, the value of D 2 was extracted from the power spectrum, allowing for 
some conhdence in the values obtained for g > — 1. Other work on known fractals indicates 
that mass-oriented approaches are superior for the analysis of low density regions, which 
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FIG. 12. Time evolution of the power spectrum of the density evaluated by averaging over 40 
successive points. At T = 20, in the interval 0,01 < A: < 1000, the slope —0,45 is obtained by 
linear regression. 

dominate when q is negative, or when 7 is positive [24]. 

For observational data, Dubrulle and Lachieze-Rey attribute the growth of the Dq{q) 
curve based on the BC or IC analysis to a lack of resolution. This argument can be considered 
here in two ways: Either the number of particles is not sufficient to properly represent low 
density areas or the analysis based on the configuration space alone is inadequate. As we 
have seen, the clusters actually form in phase, or /i, space. The projected positions on the 
line cannot account for the existence of clusters of different speeds (in other words, one 
cluster can hide another). Work on the analysis of clusters in the phase space is currently in 
progress. For the RF model with a water bag type initial condition, previous work has shown 
that Do ~ .77 [17, 31] . One approach for obtaining a better representation of low density 
regions is to use Eulerian simulations in which we follow the evolution of a distribution 
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function in phase space. Work is also in progress with this approach. 

In summary, we have seen that a one-dimensional model of an expanding universe shares 
qualitative properties with the observed universe [32], such as the bottom-up scenario, where 
due to fragmentation, a set of structures come together to form larger structures and so on. 
This hierarchical structure has a multifractal character, a property found by the results of 
our simulations. In particular, for g > 1 the generalized dimension is about 0.4 the function, 
while for g < 0 the Dq approaches about 1.4 for g about —4. This type of mutifractal behavior 
in a one-dimensional system is reminiscent of the well studied multiplicative binomial process 
[6, 24]. The function Tq exhibits two linear zones, one for g > —1, the other for g < —1. This 
suggests that the underlying geometry is bi-fractal, but caution must be excecersized because 
it reflects computations using partitions of equal size. Such a property was proposed some 
time ago by Balian and Schaeffer based on the analysis of data catalogs [33]. The work also 
directly addresses the important question of the scale on which fractal behavior is observed 
at a given epoch. We have seen that the scaling range increases with time, thus limiting the 
size of the largest objects. Beyond a horizon that grows with time, the universe as modeled 
here, is truly homogenous. Thus, in contrast with some conjectures [34-38], in the Einstein- 
de Sitter scenario the cosmological principle is verihed, but the homogeneity scale cannot 
be viewed as constant and increasingly large structures, such as the Laniakea supercluster 
[2], will develop as time progresses. 
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